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Large reversible changes of thermal conductivity are induced by mechanical stress, and the corresponding 
device is a key element for phononics applications. We show that the thermal conductivity k of ferroic 
twinned thin films can be reversibly controlled by strain. Nonequilibrium molecular dynamics simulations 
reveal that thermal conductivity decreases linearly with the number of twin boundaries perpendicular to the 
direction of heat flow. Our demonstration of large and reversible changes in thermal conductivity driven by 
strain may inspire the design of controllable thermal switches for thermal logic gates and all-solid-state 
cooling devices. 

Phonons are responsible for the transmission of sound and heat in the solid state. The ability to control 
phonon flow promises major technological developments in acoustics and thermal management 1 ' 2 , cf. 
transistors that control charge flow in electronics. Solid-state thermal switches and diodes 1 are particularly 
desirable for room-temperature cooling systems based on caloric materials 3 , as losses associated with the fluids 
that are currently employed to exchange heat between the working body and the heat sink or the heat load, 
respectively, reduce the performance of existing prototypes, and therefore cooling power 4 . Also, tunable thermal 
conductivities may be exploited for the design of advanced thermoelectric materials 5 that display high values of 
the thermoelectric figure of merit ZT = S 2 gT/k over broad temperature ranges, where S, cr, and k are the Seebeck 
coefficient, electrical conductivity, and thermal conductivity of the material, respectively, and T is the absolute 
temperature. 

In recent years, a number of strategies have been exploited to demonstrate solid-state passive thermal diodes at 
room temperature 6 " 10 , but the rectification coefficients 7c h i g h/Kiow are small. Inhomogeneously mass-loaded 
carbon and boron nitride nanotubes 6 display values of K high /Ki ow ~ 1 (1.02 for carbon nanotubes, 1.07 for boron 
nitride nanotubes), whereas few-layer graphite with asymmetric shapes display values K high /K\ ow < 2 (1.28 for 
triangular- shaped specimens 7 , 1.6 for Y-shaped specimens 8 ). Passive thermal rectification has been also demon- 
strated in bulk materials made of two oxides with different thermal conductivities 9 , but K high /K\ ow = 1.46 at best 10 . 

Tunable thermal diodes offer active control of heat flow and therefore are highly attractive for the imple- 
mentation of thermal circuits and solid-state coolers 1 . Besides various thermal diodes, control of thermal con- 
ductivity has been recently demonstrated in telescopically extended multiwall carbon nanotubes 11 (K higil /Ki ow = 
6.7), and electrostatically actuated etch-released silicon nitride microfilaments 12 (K high /K\ ow = 1.9). However, 
these intricate devices provide small areas for heat exchange, and are likely prone to fatigue during prolonged 
operation. Tunable thermal conductivity has been also demonstrated in V0 2 nanorods 13 (^hi g h/^iow = 1-6), but 
substantial Joule heating associated with the electric currents that drive the metal-insulator transition in the 
nanorods compromise rectification performance. 

Multiferroic materials constitute fertile ground to achieve sought-after large and reversible changes in thermal 
conductivity, as applied fields (electric field, magnetic field or stress) can control the number and orientation of 
phonon barriers associated with ferroic twin boundaries (TBs), whose density depends primarily on the mag- 
nitude and direction of the applied field 14 . Here we show that the thermal conductivity of ferroelastic twinned 
films can be reversibly controlled by strain to yield large thermal switching rate, due to the introduction and 
removal of twin boundaries that lay perpendicular to the direction of heat flow. Our strategy can be readily 
extended to ferroelectric and ferromagnetic twinned films whose twin boundary configuration can be controlled 
by electric and magnetic fields. Switchable domain boundaries can therefore replace previously exploited topo- 
logical defects (e.g. grain boundaries and kinks) 1516 that cannot be controlled by external fields, and lead to a new 
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paradigm for the design of thermal devices. Also, unlike competing 
mechanisms to control heat flow, we avoid the use of power- 
dissipating currents (cf. Peltier elements), and moving parts. 

We investigate coherent twin boundaries because they are most 
easily manipulated mechanically, electrically and magnetically. Our 
recent work has already shown that the morphologies and properties 
of twin patterns could be changed greatly by mechanically deforma- 
tion 17 " 19 . The important role of domain boundaries on materials 
properties 20 " 23 has been previously evidenced beyond the field of 
phononics 1 . For example, superconducting twin walls exist in W0 3 
while the matrix is non-superconducting 20 ; and ferroelastic systems 
such as Pb 3 (P0 4 ) 2 (ref. 24) and BaTi0 3 (ref. 25), and martensitic 
alloys 26 such as Ni-Ti and CuZnAl can take novel shape memory 
functionalities with the aid of domain boundary motion. Also, the 
interaction of domain wall and point defects largely determines the 
lifetime of ferroelectric memory devices 27 , while some ferroelectric 
twin boundaries are found in paraelectric bulk materials 28 , and may 
be tuned electrically. 

We simulate the effect of mechanically driven domain boundaries 
evolution on heat transfer by using non-equilibrium molecular 
dynamics (NEMD) technique 29 . The simulation details are shown 
in Method. The initial condition is a sandwich twinned structure 
containing two horizontal twin boundaries, as shown in Fig. 1. The 
surface ratio of intermediate layer to the whole sample is fixed to be 
0.7. The size of 2D simulation box is 160a X 100a in xy plane, where a 
is the lattice repetition unit. Periodic boundary condition is applied 
along the x direction and free boundary condition is used in the y 
direction. We apply different kinds of strain state to examine the 
twinning evolution in the mechanically driven system. The strain 
tensor is composed of [s xx , s yy , y xy ] T , where s xx is tensile strain along 
x direction, s yy is tensile strain along y direction, and y xy is shear 
strain. In our work, the deformation is performed at different s yy / 
y xy ratio with s xx = 0. The simulations details are shown in Methods 
and Supplementary Section 1. 

We first apply a simple shear deformation to the sample with 
strain tensor as [0, 0, y xy ] T . Figure 2a shows the response of shear 
stress with shear strain (y xy ) in a cycle. Under loading (black curve) 
the sample yields when y xy reached —0.6%. The sample was unloaded 
to zero-strained state from y xy = 1.6% (red curve). Figure 2b shows 




Figure 1 | Thin-film model. A two-dimensional sandwich model with two 
fixed horizontal twin boundaries (HTBs). The lattice unit is in the shape of 
a parallelogram with tilt angle of 4 degrees. The middle layer has an area 
ratio of 0.7 to the whole sample. The strain tensor in 2D system can be 
described as [e^, s yp y xy ] T . In our work, the deformation was performed at 
different s yy /y X y ratio with = 0. Heat flow will be applied along x 
direction. Atoms are colored according to values of their centro- symmetry 
parameters in range of 0 ~ 1 (ref. 41). 



the variation of k during the loading/unloading loop. Here k is given 
in relative unit, r.u., related to the model potential we used here. We 
have not calibrated the absolute values to those which are obtained 
experimentally. In elastic regime, the magnitude of k does not change 
significantly (—140 to 150 r.u.) under stress. When loading into the 
plastic regime, k undergoes an abrupt drop to —90 r.u., almost one 
half of the initial value. This magnitude is sustained under further 
stress. When the system is unloaded, k increases again and shows 
hysteretic behavior during the loading/unloading cycle. Our work 
here shows that by mechanically controlling the microstructures 
rather than the intrinsic properties of the single domain state, it is 
possible to generate two logic states, one with high heat conduction 
(such as points c, f in Fig. 2b) and the other with low heat conduction 
(such as points d, e in Fig. 2b), which can be used for thermal 
information storage and thermal switching 13 ' 30 . 

We examined the evolution of twin pattern for several typical 
strains (marked as (c), (d), (e), (f) in Fig. 2b) upon loading/unloading, 
as shown in Figs. 2c-f. When the system deforms plastically, new 
horizontal twin layers are induced by the applied shear strain, accom- 
panied with the formation of a certain amount of vertical twin 
boundaries (VTBs). VTBs nucleate from one horizontal twin bound- 
ary (HTB), propagate and terminate in another horizontal twin 
boundary. The new-formed VTBs and HTBs superimpose and 
finally evolve into a much more complicated twin pattern. This pat- 
tern then shows much reduced thermal conductivity. Moreover, the 
structure of twin pattern is quite stable and the twinning morphology 
will not undergo large changes even under the applied temperature 
gradient in NEMD calculations (see Supplementary Section 2). Upon 
unloading, horizontal twins are preserved and produce a permanent 
deformation of the sample. Vertical twins, on the other hand, are 
unstable and vanish gradually when removing the external load. This 
leads to an increase of k. 

These results show that the reduction of thermal conductivity 
results from the existence of vertical twins with an orientation per- 
pendicular to the heat flow. To further probe the exact role that VTBs 
and HTBs play in thermal transfer, we examined the variation of both 
VTBs density (Pvtb) and HTBs density (Phtb) in loading/unloading 
cycle. Twin densities in our 2D system are calculated as the ratio of 
twin boundary area to total volume of system, which is in unit of a" 1 . 
Figure 2g shows that the production and annihilation of vertical 
twins are correlated with changes of the thermal conductivity, i.e. 
the formation of VTBs reduces thermal transport. As the stress 
reaches the upper yield point (y xy = 0.8%), VTBs begin to nucleate, 
which corresponds to an abrupt drop of thermal conductivity. The 
magnitude of Pvtb reaches to the maximum value at lower yield 
point (y xy = 1.0%) and is kept around that level under further shear- 
ing. Correspondingly, the thermal conductivity is also almost 
unchanged in this strain range of 1.0%-2.0%. In contrast, the pres- 
ence of horizontal twins does not influence the heat transport, and 
thermal conductivity is independent of the density of HTBs (Fig. 2h). 
To distinguish the two twin structures with different heat transfer 
properties, we term the system containing only HTBs as twinning 
pattern 1 (TP1) and the system containing VTBs (purely vertical twin 
boundaries or mixed twin structure) as twinning pattern 2 (TP2). For 
the present case, the magnitude of k in TP2 is halved compared with 
TP1. The above calculations agree well with the experimental obser- 
vation that the existence of twins in ferroelastic Gd 2 (Mo0 4 ) 3 reduces 
thermal conductivity largely in comparison with a monodomain 
sample 31 . 

We have shown that VTBs density is the dominant factor for heat 
transportation properties. Since the twin patterns are sensitive to 
external strain state, it is possible to manipulate VTBs densities by 
applying different strain fields to tune thermal conductivity. When 
replacing the simple shear to a mixture of loading state [0, s yy , y xy ] T , 
i.e. applying certain amount of tensile strain (s yy ) along y direction 
besides shear strain (y xy ), a much higher magnitude of resolved shear 
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Figure 2 | Strain-control of thermal conductivity in twinned films, (a) Stress-strain curve in a simple shear loop, (b) Variation of thermal 
conductivity with shear strain, (c)-(e) Representative atomic images corresponding to the different strain-scenarios marked in (b). Twin pattern 1 (TP1) 
in absence of vertical twin boundaries ( VTBs) has a lower thermal conductivity than twin pattern 2 (TP2) with VTBs. Density of (g) VTBs and (h) HTBs as 
a function of shear strain, in units of a~\ Atoms are colored according to values of their centro- symmetry parameters 41 . 
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Figure 3 | Enhancing thermal-conductivity tunability using multiple strains, (a) Thermal conductivity variation with shear strain under different e yy ly xy 
ratio, (b) Variation of VTBs density with shear strain, (c)-(e) Representative atomic images corresponding to the different strain scenarios marked in (b). 
Atoms are colored according to values of their centro- symmetry parameters 41 . 



stress will be imposed on the VTBs planes and the density of VTBs 
can be enhanced greatly. Figure 3a and 3b show the variation of 
thermal conductivity and VTBs density with applied shear strain 
under different s yy /y xy ratio, respectively. The VTBs density increases 
continuously when a tensile strain is applied. Correspondingly, the 
thermal conductivity is reduced dramatically. Figures 3c-e show the 
atomic configurations at y xy = 1.8% under different s yy /y xy ratios. We 
find that the application of tensile strain effectively suppresses the 
formation of HTBs, and enhances the formation of VTBs. 
Furthermore, thermal conductivity can be reduced to 46 r.u. (see 
orange point in Fig. 3a), which is more than 30 smaller than the bulk 
values 1488 r.u. (see Supplementary Section 3). 

The atomic mechanism of the reduction of thermal conductivity 
by VTBs is seen from the analysis of phonon spectra and vibrational 
eigenmodes 32 . We choose a sample with the system size of 80a X 50a 
to calculate the phonon density of states (DOS), keeping all other 
simulation conditions the same. Figure 4a shows the DOS for differ- 
ent values of Pvtb extending to 27 THz. With increasing p V TB the 
spectral weight at high frequencies (15-27 THz) decreases. This 



implies the twin boundary-phonon scattering by VTBs are at high 
frequencies while low frequency phonons are much less affected by 
VTBs (increasing is due to the normalization of the phonon spectra). 
To quantitatively describe the impact of VTBs on the phonon trans- 
port, we calculate the participation ratio px of each mode, which is 
defined as 33 ' 34 

i a 

where N is total number of atoms and s i0Ly x is the ath eigenvector 
component of each eigenmode X. The participation ratio describes 
the fraction of atoms participating in a mode and provides the 
information on localization effect of each mode. The value of px 
varies from the 0(1/ N) for localization to 0(1) for derealization. 
We find that more phonons become localized when p VTB increases 
and correspondingly thermal conductivity is reduced (Fig. 4b). In 
particular, significant phonon localization occurs at high frequencies 
(around 18-27 THz). This is related to the highly distorted atomic 
structure at VTBs. During the twin patterns evolution, atoms inside 
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Figure 4 | Phonon analysis of the influence of VTBs on thermal 
conductivity. Variation of (a) phonon spectrum and (b) participation 
ratio of systems under different p VTB during shear plasticity, (c) Bulk 
phonon dispersion relation. Red and black lines refer to the longitudinal 
and transverse branches, respectively. 

twin boundaries displace and lead to phonon scattering at the smal- 
lest length scale of k = 0. 5 — In/ a and corresponding frequencies 
around 18-27 THz (see phonon dispersion relation in Fig. 4c). In 
addition, Fig. 4c clearly shows that phonons of 18-27 THz, beyond 
the transverse mode, are specific to the longitudinal mode. Therefore, 
we see clearly that the emergence of VTBs strongly scatters the trans- 
port of the longitudinal phonon transport that in turn severely sup- 
presses the longitudinal thermal transport. 

In kinetics theory of phonon transport, the thermal conductivity is 
proportional to the mean free path of phonons as k = ACv/3, where 
A is the mean free path, C is the heat capacity and v is the acoustic 
velocity. In our system, the mean free path is determined by two 



kinds of phonon scattering events, involving (a) phonon-phonon 
collision and (b) interactions of phonons and twin boundaries. The 
effective A is 



A Ap_p Ap_tb ' 

where A P . P and A P . TB are the characteristic lengths of the bulk pho- 
non-phonon scattering and the phonon-TBs scattering. The mag- 
nitude of Ap.p in ferroelastic materials is usually on the order of 0.1 
micrometers at low temperature 31 , which is much larger than A P . TB 
here. Thus, A P . TB plays the dominant role in determining the total 
mean free path A. For a constant A P . P , 1/k should be proportional to 
l/Ap_ TB (1/k oc 1/A P _ TB ). On the other hand, here A P . TB should be 
equal to the average twin boundary spacing (d), which is inversely 
proportional to the density of VTBs d = 1/Pvtb- Therefore, we can 
finally obtain an approximate relationship of 1/k and Pvtb as 

^OCPVTB, (3) 

To test this proportionality, we extract data within the strain range of 
1.2-1.8% for different strain states in Fig. 3a and plot them in Fig. 5a. 
The data fit the linear relationship between 1/k and Pvtb in Eq. 3 very 
well. It also manifests that the correspondence of 1/k and p V TB is not 
affected by the strain state we applied. However, it should be noticed 
that the intersection point of linear line with)/ axis (1/k 0 ) does not 
refer to the bulk value (l/?c bu i k ) because systems with shorter box 
lengths (L x ) will underestimate the value of 7c bu i k due to the effect of 
periodic boundary condition on mean free path of phonons 
(Supplementary Fig. S3). The inset in Fig. 5b shows the dependence 
of simulation length on 1/k 0 and the value converges to l//c bu i k ( = 6.7 
X 10" 4 r.u.). Moreover, the linear relationship described in Eq. 3 can 
be further verified in the twinned system with larger lengths. The 
slope 1/k to Pvtb increases correspondingly and reaches 0.14 (orange 
line) in our calculations. Referring the slope to the bulk system as 
reference, this slope is predicted to become even steeper when 
removing the limitation of small length scales in our present MD 
simulations. Our results are compiled in Fig. 5b where the reference 
point of the untwinned materials is 1/ftbuik and all other values relate 
to nanostructured materials with a defined number of vertical twins 
during cold shearing and different sample sizes. Twin densities of 
0.1a" 1 (i.e. twin walls separated by 10 unit cells) are plausible in view 
of wall separations of 15 nm (some 40 unit cells) in SrTi0 3 when the 
walls are thermally induced 35 . Cold shearing leads to much higher 
wall densities so that device applications become very likely in the 
light blue triangle in Fig. 5b. At a likely density of 0.1a" 1 in light blue 
triangle, 1/k = 0.015 r.u. which is larger by a factor of 22 than the 
equivalent bulk value. An upper limit of the effect can be estimated 
when the twin density increases to 0.16a" 1 . The simulate contrast of 
the thermal conductivity reaches a factor K: h i g h/^buik = 30 in this case. 
However, it should be emphasized that the thermal switching rate is 
fixed for a given sized system. For example, K high /Ki ow ~ 2.9 when L x 
= 160a (see the black open symbols in Fig. 5b). Based on these 
results, we are able to approach the goal of tuning thermal properties 
in the framework of "domain boundary engineering" and estimate a 
maximum possible reduction of thermal conductivity by applied 
stress where the density of VTBs can be manipulated very precisely 
by external strain (or stress). This task becomes even easier with 
recently developed nanotechnology which fabricates materials with 
nanotwinned structure 36 . Since VTBs reduce thermal transport, we 
can hence tune thermal conductivity via the applied strain to the 
sample, which, in turn, can be manipulated by electric or magnetic 
fields in multiferroic materials. 

Overall, our work shows how to achieve large and reversible 
changes in the thermal conductivity of multiferroic twinned films 
using applied fields. By controlling magnitude and direction of the 
applied field, we manipulate the density of twin boundaries that lay 
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perpendicular to the direction of heat flow and act as phonon bar- 
riers. Our strategy may be exploited for the implementation of fast 
controllable thermal devices, where heat- exchange surfaces can be 
largely increased by using multilayer geometries. 

Methods 

Our potential is constructed to be generic for all ferroelastic phase transitions by 
choosing the shear angle as "order parameter" and can be used to study twinning and 
mobility of twin boundaries 14 in ferroelastic materials. Comparing this configuration 
with the ferroelastic transition from cubic to tetragonal in SrTi0 3 at 105 K, we focus 
on the modeling of a plane perpendicular to the twin plane (i.e. our axes are [101] and 
[1 01] in the Pm 3 m setting). Equivalent planes can be found in all ferroelastic 
systems. The Landau potential of such materials is implemented in the model to 
represent the interactions in a mono-atomic, two-dimensional lattice. The potential 
energy U(r) contains three parts: first-nearest atomic interactions of 20(r — l) 2 , 
second-nearest interactions — 10(r — a/2) 2 + 2000(r — \fl) A , and third-nearest 
interactions — (r — 2) 4 , where r is the atomic distance. Details of this potential are 
described in our previous work 17 " 19,21 . The equilibrium unit cell is in the shape of a 
parallelogram with shear angle of 4 degrees from square. 



To study the structure evolution upon the external load in the sandwich sample, the 
top and bottom three atomic layers were fixed rigidly as the loading grip. We apply 
different kinds of strain state to examine the twinning evolution in the mechanically 
driven system. The strain tensor is composed of [e xx , s yy , y xy Y, where e xx is tensile 
strain along x direction, & yy is tensile strain along y direction, and y xy is shear strain. In 
our work, the deformation is performed at different s yy /y xy ratio with s xx = 0. The 
dynamic loading is taken at T ~ 0.2 T m by using Nose-Hoover thermostat 3738 , where 
T m is melting point. Under such low temperature, any atomic diffusion process can be 
greatly suppressed. All the calculations are performed using the LAMMPS code 39 . The 
atomic configurations are displayed by AtomEye 40 and the colors are shown 
according to the centro- symmetry parameters 41 . 

For each loading step, we used the NEMD technique to calculate thermal con- 
ductivity in a given configuration 29 . The idea is to apply a heat flux in the system along 
the x direction, which will result in a temperature gradient. When the heat transfer 
becomes a steady flow, the thermal conductivity k along x direction is calculated via 
Fourier law as k = —J x /(dT/dx), where J x is the heat flux from heat source to heat sink 
and T is temperature. The induced temperature gradient at steady- flow ranges from 
0.28 T m to 0. 1 3 T m , corresponding to the heat source and the heat sink (Supplementary 
Fig. SI). Note that real films, with larger unit cells and multiple atomic species, will 
exhibit much higher melting points and therefore permit room temperature 
operation. 
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Finally, we note that the present simulations were performed in 2D systems. 
Although our recent work shown that 2D simulations can reproduce the same 
thermodynamics properties in ferroelastic systems as in 3D simulations 42 , a study of 
three-dimensional materials would be highly desirable but also extremely costly in 
computer time. 
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